# Analysis and figures for:
# Unions, Class Identification and Policy Attitudes
# Date: 03/31/22

# R version 4.1.2


### Change path to directory where replication files are stored. ###
setwd("~/Dropbox/Class/SubjectiveClass/JOP_Replication/")


library(haven)
library(ggplot2)
library(grid)
library(gridExtra)
library(dplyr)
library(forcats)
library(ggstance)


# Set theme options.
color.background <- "#F0F0F0"
color.grid.major <- "#DBDBDB"
fte_theme <- function() {
  theme(plot.background = element_rect(fill = color.background)) +
    theme(panel.background = element_rect(
      fill = color.background,
      colour = color.background)
      ) +
    theme(panel.grid.major=element_line(
      color = color.grid.major,
      size = .25)
      ) +
    theme(panel.grid.minor = element_blank()) +
    theme(axis.ticks = element_blank()) +
    theme(axis.title.x = element_text(vjust = -0.5), 
          axis.title.y = element_text(vjust = 1.2)) +
    theme(plot.title = element_text(hjust = 0.5)) +
    theme(legend.text = element_text(size = 10)) +
    theme(legend.title = element_blank()) +
    theme(legend.background = element_rect(
      fill = color.background, 
      color = "black", 
      size = 0.1)
      ) +
    theme(legend.key = element_rect(color = color.background)) 
}



### FIGURE 1. ###
# Interaction effects.
# Load data.
# The 'rename_all' code removes underscores from variable names.
pr.inter.ML <- read_dta("pr_objXuni_ML.dta") %>% 
  rename_all(~ gsub("_", "", .))

# Create factor variables.
pr.inter.ML <- pr.inter.ML %>%
  # Outcome.
  mutate(outcome = case_when(predict == 1 ~ "Subj. lower class",
                             predict == 2 ~ "Subj. working class",
                             predict == 3 ~ "Subj. mid. class",
                             predict == 4 ~ "Subj. upper class")) %>%
  mutate(outcome = fct_reorder(outcome, predict)) %>%
  # Objective class.
  mutate(objcls = case_when(m1 == 1 ~ "Obj. lower",
                            m1 == 2 ~ "Obj. working",
                            m1 == 3 ~ "Obj. lower mid.",
                            m1 == 4 ~ "Obj. upper mid.", 
                            m1 == 5 ~ "Obj. upper")) %>%
  mutate(objcls = fct_reorder(objcls, m1)) %>%
  # Union mem.
  mutate(union = case_when(m2 == 0 ~ "Non union",
                           m2 == 1 ~ "Union")) %>%
  mutate(union = fct_reorder(union, m2))


fig.inter.ML <- ggplot(pr.inter.ML, 
                       aes(y = objcls, x = margin, 
                           color = union, shape = union)) +
  geom_errorbarh(aes(xmin = cilb, xmax = ciub), 
                 size = 1, color = "gray50", width=0) +
  geom_point(size = 2) +
  labs(x = "Probability of subjective class ID", 
       y = "") +
  facet_wrap( ~ outcome, ncol = 2) +
  fte_theme() + 
  theme(legend.position = "bottom") +
  scale_shape_manual(values = c(1, 2))

ggsave("Figure_1.pdf", width=6.5, height=5.0)



### FIGURE 2. ###
# Political attitudes.
# Load data.
scls_coef2 <- read_dta("scls_coef2.dta") %>%
  filter(!(parm %in% c("1b.region4", "2.region4", 
                       "3.region4", "4.region4")))

# Reorder factors.
scls_coef2$var <- factor(scls_coef2$var, 
                         levels = c("Subj. lower",
                                    "Subj. working (ref.)",
                                    "Subj. middle",
                                    "Subj. upper",
                                    "Obj. lower",
                                    "Obj. working (ref.)",
                                    "Obj. lower mid.",
                                    "Obj. upper mid.",
                                    "Obj. upper",
                                    "Party ID",
                                    "Ideology",
                                    "Union mem.",
                                    "Female",
                                    "White (ref.)",
                                    "Black",
                                    "Other"))
scls_coef2$var <- fct_rev(scls_coef2$var)

scls_coef2$depvar <- factor(scls_coef2$depvar, 
                            levels = c("Gov. reduce ineq.",
                                       "Gov. help poor",
                                       "Gov. help sick",
                                       "Gov. do more"
                            ))


# Plot.
fig.subjpol2 <- ggplot(scls_coef2, 
                       aes(y = var, x = b)) +
  geom_vline(xintercept = 0, color = "gray", 
             linetype = "dashed") +
  geom_errorbarh(aes(xmin = ll, xmax = ul), 
                 size = 1, color = "gray50", width=0) +
  geom_point(size = 1.5, color = "seagreen3", shape = 1) +
  labs(x = "", 
       y = "") +
  facet_wrap( ~ depvar, nrow = 1) +
  fte_theme() +
  theme(plot.background = element_rect(fill = "white")) 


ggsave("Figure_2.pdf", width=7.0, height=4.0)




### Appendix figures. ###

# Load data for figures A1 and A2.
gssdat <- read_dta("GSS_SubClass_clean.dta")
gss.0616 <- filter(gssdat, year>=2006 & year<=2016, 
                   !is.na(egp5))

### FIGURE A1. ###
# Bar plot of income and education by EGP class categories.
# Labels.
egp5_labs <- c("1"="Upper Class", "2"="Upper Middle Class", "3"="Middle-Class Serv.", 
               "4"="Middle-Class Man.", "5"="Working Class")

# Income and education within EGP class categories.
figEGP.inc <- ggplot(filter(gss.0616, !is.na(egp5), !is.na(inc5)), 
                     aes(factor(inc5), weight = wtss)) +
  geom_histogram(aes(y = ..count.. / tapply(..count..,..PANEL.., sum)[..PANEL..]), 
                 stat = "count", fill = "#66CC99") +
  labs(y = "", x = "Respondent Income") +
  facet_wrap( ~ egp5, ncol = 2, labeller = labeller(egp5 = egp5_labs)) +
  fte_theme() 


figEGP.edu <- ggplot(filter(gss.0616, !is.na(egp5), !is.na(edu5)), 
                     aes(factor(edu5), weight = wtss)) +
  geom_histogram(aes(y = ..count.. / tapply(..count..,..PANEL.., sum)[..PANEL..]), 
                 stat = "count", fill = "#619cff") +
  labs(y = "", x = "Respondent Education") +
  facet_wrap( ~ egp5, ncol = 2, labeller = labeller(egp5 = egp5_labs)) +
  fte_theme() 

# Combine plots.
pdf("Figure_A1.pdf", width=6, height=3)
grid.arrange(figEGP.inc, figEGP.edu, ncol=2)
dev.off()



### FIGURE A2. ###
# Histograms of latent objective class variable.

# Get min and max values.
oc.min1 <- aggregate(objcls2 ~ objcls5, gssdat, min)
oc.max1 <- aggregate(objcls2 ~ objcls5, gssdat, max)

# Plots.
objhist.gss <- ggplot(gssdat, aes(x = objcls2)) + 
  geom_vline(xintercept = c(oc.min1$objcls2, oc.max1$objcls2), 
             color = "black", lty = "dashed") + 
  geom_histogram(binwidth = 0.08, fill = "#f8766d") +
  labs(x = "Latent objective class") +
  fte_theme()

ggsave("Figure_A2.pdf", width=6.0, height=5.0)



### FIGURE A3. ###
# Plot predicted probabilities.
# Load data.
# The 'rename_all' code removes underscores from variable names.
pr.scls.gss <- read_dta("pr_scls.dta") %>% 
  rename_all(~ gsub("_", "", .))

# Reorder factors.
pr.scls.gss$vlab <- factor(pr.scls.gss$vlab, 
                           levels = c("Obj. lower",
                                      "Obj. working",
                                      "Obj. lower mid.",
                                      "Obj. upper mid.", 
                                      "Obj. upper",
                                      "Non union",
                                      "Union"))

# New category by variable.
pr.scls.gss <- pr.scls.gss %>%
  mutate(vcat = case_when(var >= 1 & var <= 5 ~ "Objective class",
                          var >= 6 & var <= 7 ~ "Union mem."))


# Labels for outcomes.
pr.scls.gss <- pr.scls.gss %>%
  mutate(outcome = case_when(predict == 1 ~ "Subj. lower class",
                             predict == 2 ~ "Subj. working class",
                             predict == 3 ~ "Subj. mid. class",
                             predict == 4 ~ "Subj. upper class")) %>%
  mutate(outcome = fct_reorder(outcome, predict))


# Plot.
fig.scls.gss <- ggplot(pr.scls.gss, aes(y = factor(var), x = margin)) +
  geom_vline(xintercept = 0, color = "gray", 
             linetype = "dashed") +
  geom_hline(yintercept = 5.5, color = "white", size = 2) +
  geom_errorbarh(aes(xmin = cilb, xmax = ciub), 
                 size = 1, width=0) +
  geom_point(size = 1.5, color = "#619cff", shape = 1) +
  facet_wrap( ~ outcome, ncol = 2) +
  labs(x = "Predicted probability of subjective class ID", 
       y = "") +
  fte_theme() +
  theme(plot.background = element_rect(fill = "white")) +
  scale_y_discrete(labels=c("1" = "Obj. lower",
                            "2" = "Obj. working",
                            "3" = "Obj. lower mid.",
                            "4" = "Obj. upper mid.", 
                            "5" = "Obj. upper",
                            "6" = "Non union",
                            "7" = "Union"))


ggsave("Figure_A3.pdf", width=6.0, height=4.5)



### FIGURE A4. ###
pr.inter.WC <- read_dta("pr_objXuni_WC.dta") %>% 
  rename_all(~ gsub("_", "", .)) 

# Create factor variables.
pr.inter.WC <- pr.inter.WC %>%
  # Objective class.
  mutate(objcls = case_when(m1 == 1 ~ "Obj. lower",
                            m1 == 2 ~ "Obj. working",
                            m1 == 3 ~ "Obj. lower mid.",
                            m1 == 4 ~ "Obj. upper mid.", 
                            m1 == 5 ~ "Obj. upper")) %>%
  mutate(objcls = fct_reorder(objcls, m1)) %>%
  # Union mem.
  mutate(union = case_when(m2 == 0 ~ "Non union",
                           m2 == 1 ~ "Union")) %>%
  mutate(union = fct_reorder(union, m2))


# Plot.
fig.inter.WC <- ggplot(pr.inter.WC, 
                       aes(y = objcls, x = margin, 
                           color = union, shape = union)) +
  geom_errorbarh(aes(xmin = cilb, xmax = ciub), 
                 size = 1, color = "gray50", width=0) +
  geom_point(size = 2) +
  labs(x = "Probability of subjective working class ID", 
       y = "") +
  fte_theme() + 
  theme(legend.position = "bottom") +
  scale_shape_manual(values = c(1, 2))


ggsave("Figure_A4.pdf", width=6.5, height=5.0)


